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Introduction 


In this study we shall discuss ways of accelerating the approach to a 
steady state for the Euler equations* We first consider the time-dependent 
equations* However, since we are only interested in the steady state we shall 
feel free to alter the equations in any way that does not affect the steady 
state* Since the Euler equations constitute a nonlinear system of equations, 
it usually can not be proven that the modified system even approaches a steady 
state* Instead, most of the analysis will refer to linearized versions of the 
equations* Assuming that the steady state is unique once the modified 
equations reach a steady state, it must be the steady state of the original 
equations* We shall concentrate on the inviscid equations, though most of the 
techniques are also applicable to the Navier-Stokes equations* 

We assumed that a body— fitted curvilinear grid has been constructed from 
some package. All that we require is the (x,y) coordinates of each zone* 
The finite difference equations will be derived using a finite volume 
approach* Though the cells can be of any shape in such an approach we shall 
assume that all cells are quadrilaterals* Near the trailing edge some zones 
may degenerate and the finite volume approach still holds* The use of mapping 
formulas rather than finite volumes yields essentially the same finite 
difference formulas* The finite volume approach naturally associates the 
value of the dependent variables with zonal averages* These averages are 
associated with values at the cell center* Hence, all boundaries are placed 
along cell faces* With a mapping formulation the natural implementation would 
be to place the variables at grid nodes* Furthermore, using a finite volume 
approach, the Jacobian of the transformation is identified with the area of a 
cell* This leads to a slightly different formulation for the Jacobian than 
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would be natural for a finite difference approach. The finite volume 
formulation leads to a straightforward generalization for axisymmetric 
coordinates. 

For convenience of notation we shall only consider two-dimensional flows. 
The Euler equations can then be written as 

Wj. + f + gy = 0, (1.1) 


T 

w = (p, pu, pv,E) , 

2 T 

f = (pu, pu + p, puv, puh) , 

( 1 . 2 ) 

2 T 

g = (pv, puv, pv + p, pvh) , 



Integrating (1.1) inside a cell and using the divergence theorem we get 


|— // wdV + / (fdy - gdx) = 0, (1.3) 

where we have chosen D to be a quadrilateral. Let two adjacent sides of the 
quadrilateral be given by 5 = constant and n = constant. The component of 
velocity perpendicular to the curve ? = constant, denoted by q, is given by 
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while the velocity component along the curve 5 
given by 

_ •_!] n 


constant, denoted is 

(l.Ab) 


Similarly, letting r be the velocity component perpendicular to q = 
constant and rj^j^ parallel to n = constant, we have 


r = Hx u + riy V = (-y^ u + v)/J, 


(1.4c) 


Til = (y^ V + u)/J, 


(1.4d) 


r is proportional to qn» only if the grid is orthogonal, i.e. , 
x- y + y- X =0. 

Letting w^ denote the cell average 


''a “ 


/ / wdv 
// dv 


(1.5) 


we replace (1.3) with 


(\j •y) 
dt '■ A 


4 

+ I 

i=l 


(fdy - gdx) = 0, 


( 1 . 6 ) 


where are the sides of the quadrilateral and V is its area. We 

evaluate the line integrals using the midpoint rule. Normalizing the mesh so 
that the length of each quadrilateral is one, we replace (1.6) by the 


approximation 
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dt 


(w -V) = - I (fy 
5=const 



+ I [iyp - gx ). 

n=const 


(1.7) 


With second-order accuracy we can find w at the cell face by averaging w 
from the cell centers* This averaging is done in FL052ST without accounting 
for the difference in volumes of neighboring cells. Given w at the cell 
face, one can then calculate the fluxes* It is preferable to average w 
rather than average the fluxes to help couple the even and odd points* For 
efficiency the pressure is also averaged rather than computed from w* The 
error in this procedure is of the same order as the error in the scheme* 
Several numerical checks have confirmed that this pressure averaging does not 
introduce any additional errors* 

In the following sections we discuss ways to integrate (1*8) in time and 
to accelerate the convergence to a steady state* To measure this acceleration 
we need a way of deciding when a numerical steady state has been reached* In 
[8] the measure used was Ap/At* However, once acceleration techniques are 
used, it is important to measure the true residual, i*e*, the right-hand-side 
of (1*8)* For exterior problems the mesh is finer near the body and is 
coarser in the far field* As such the maximum residual usually occurs in the 
far field where the volumes are large* This is true even though the true flux 
is small in the far field even relative to the zonal areas* An alternative is 
to normalize the true residual by the cell area which tends to emphasize the 
zones near the airfoil* 
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2. Runge-Kutta Hethods 

In 1960 Lax and Wendroff [11] introduced a one-step method which was 
second-order in space and time. Richtmyer [15] suggested a two-step method 
which is linearly equivalent to the Lax-Wendroff scheme. The advantage of the 
two-step scheme is that only flux evaluations are necessary rather than matrix 
multiplications . 

A disadvantage of these two-step schemes is that the two steps are 

different, complicating the coding. A more important disadvantage is that the 
Courant number is one, even though two steps are used. The leapfrog method 
requires only one step and has the same stability requirement. Hence, the 

leapfrog scheme is twice as efficient as a two-step Lax-Wendroff method. 
However, the leapfrog scheme is nondissipative and hence not useful for 

shocked flows. Nevertheless, we are interested in multistage schemes which 
are more efficient than the standard two-stage schemes. Graves [7] used a 
three-stage Runge-Kutta scheme for the Navier-Stokes equations. Van der 
Houwen [21], [22] analyzed multistage schemes for both hyperbolic and 

parabolic problems. A fourth-order Runge-Kutta scheme for the Euler equations 
was popularized in [8]. For simplicity we only consider two-space dimensions 
although the extension to higher dimensions is straightforward. Let 

w = f + g . (2.1) 

t X ®y 

An N stage Runge-Kutta scheme is given by 


w^^^ = w^®^ + f 


( 2 . 2 ) 
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with = w*', , Ojj = 1. We note that this method is, in 

general, only first-order accurate in time even when the amplification matrix 
agrees with that of a higher-order method. In many cases one can achieve 
second-order accuracy in time. When one is only Interested in the steady 
state solution, this is not a drawback. For time-dependent problems one may 
wish to have higher accuracy. One can then replace (2.2) with a true higher 

order Runge-Kutta method. The advantage of (2.2) is that only two levels of 

storage are required, i.e. , w^^^ and w^^^ at each stage. 

A fundamental question is how to choose the parameters Oj^, k=l,»*»,N-l. 

To simplify the question, we only consider the linear problem with constant 

coefficients. We then Fourier transform (2.2) and consider the following 
amplification matrix 


G=1+BjZ+P2Z + 




(2.3) 


is the Fourier transform of At(D w + D w). 

X y 

At 


Hence, if we use central 


differencing in x and y, then z = i • • X(A sin 9 + B sin (|)) where 

X(A) denotes an eigenvalue of A. 0, <|> are the Fourier variables. Ax = Ay, 
A = , B = . The parameters 0^ are given by 


1 , 


^2 " °N-1’ 


(2.4) 


^k “ ^k-1 ®N-k+l* 


k = 3,‘«‘,N. 



-7- 


One possibility is to choose the 3^ so that At is as large as possible. 

If the iteration is close, in some sense, to being time accurate, this implies 

that we advance in time as fast as possible. This requires that |z| be as 

* 

large as possible while requiring that G G _< 1. In general, this At is 
only required for some frequencies which depend on the scheme and the matrices 
A and B. Using central differences the values of z lie along the 
imaginary axis while for an upwlnded scheme, z will lie on some curve in the 
complex plane. 

Using a central difference scheme the maximum time step is given by 


Ax — max p(A sin 0 + B sin (f) ’ 

0 ,<|, 

where p(A) is the spectral radius of A and K is the maximum that z can 
achieve in (2.3). Vichnevetsky [25] (see also [16]) has shown that when z 
is purely imaginary that 

|z| ^ N-1 for an N stage scheme. (2.6) 


This maximum can be achieved and the formula is given by [16], [21] 


G = Pjj(z) = i‘ 




For N odd, this formula coincides with that given by Van der Houwen [21] and 
the formula is second-order accurate. For N even, the formula is only 
first-order accurate. 
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In practice, 
(see section 4). 


an artificial viscosity is needed to stabilize the scheme 
This can be modelled, in one dimension, by 


w 


t 


w 


X 


- ew 

xxxx 


( 2 . 8 ) 


Using central differences, z now corresponds to z = 1 sin 0 - 16 e sin^ -|- . 
|z| now lies in the left half plane. Plotting the stability region of (2.3), 
(2.7) we see that the stability region contains the imaginary axis up to N-1 
and also includes portions of the right-half plane. The stability region also 
includes part of the negative real axis up to z of about two to three 
(depending on N). Hence, as long as e is sufficiently small the previous 
results are valid. At 0 = it , z = — 16 e and so the stability is governed only 
by the artificial viscosity. A similar analysis holds when the Navler-Stokes 
terms are added to the differential equation. 

In the previous analysis we assumed that the artificial viscosity is 
reevaluated at each stage. In practice the computation of the artificial 
viscosity is expensive and so the same artificial viscosity is used for each 
stage within a cycle. Hence, (2.8) is approximated by 

4t[D„ d! „«»1. (2.9) 

Let z be the Fourier transform of At • w and let y be the Fourier 
transform of At (artificial viscosity). Then the amplification factor is 


G = 1 - (z+y)[8j + z + ••• + z^"^] + 3^^ z** + y. 


( 2 . 10 ) 
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Therefore, the previous analysis does not hold since G depends on both z 
and y and not on one complex number. In various computational trials the 
optimal parameters (2.7) worked when the artificial viscosity was reevaluated 
at each stage. However, for N greater than four, the optimal parameters 
were not stable when the viscosity was frozen. It is not clear how much of 
these details are dependent on the exact formulation of the artificial 
viscosity. 

When the artificial viscosity is reevaluated at each stage, then the use 
of a N stage formula used twice is equivalent to some 2N stage formula. 
Hence, it is always more optimal to use a higher stage formula. However, when 
the artificial viscosity is frozen at each stage, this no longer need be true. 
With the second-order central difference scheme we have seen, (2.6), that the 
Courant number for an N stage scheme has a maximum of N-1. Thus, the 
efficiency per stage is (N-l)/N. Hence, for N large, we approach the 
efficiency of the leapfrog method. At N = 10, we already have 90% efficiency 
and there is not much benefit in using a higher stage method. However, since 
this scheme requires the reevaluation of the artificial viscosity at each 
stage, these higher stage schemes are not efficient. An alternative is to 
evaluate the viscosity M times within an N stage scheme. The relationship 
between M and N to achieve maximum stability is not known. Furthermore, 
in many cases, robustness, i.e., including sections of the left-half of the 
complex plane, is more important than maximal time steps. 

Until now we have assumed that the stability condition is mainly governed 
by z near the imaginary axis, i.e., central differences with a small 
artificial viscosity. Equations with variable coefficients or nonlinearities 
require us to consider perturbations off the imaginary axis even without the 
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presence of viscous effects, see [15]. For parabolic equations the 
appropriate z are along the negative real axis. The optimal parameters for 
this case were first considered in [17]. This is given by 

G = T^(l + z/N^), 0 < z < -2N^ (2.11) 

for an N stage scheme. However, this formula is not stable for large 

negative z if z is slightly off the real axis. These perturbations can be 
caused by variable coefficients, lower order terms or boundary conditions. 
One type of perturbation is considered in [22]. Another perturbation of 

(2.11) is 

G = (1 - e)Tjj(y) + e y Tj^(y)/N^, with y = 1 + z/N^. (2.12) 

The stability range of (2.12) is slightly smaller than but contains a large 

region off the real axis when e > 0. In both cases the stability region for 
an N stage scheme is proportional to N^. Hence, the more stages that are 
used the more efficient the method is. This is in contradistinction to the 
hyperbolic case where the efficiency per stage quickly approaches its upper 
bound which is given by the leapfrog scheme. Hence, when the number of stages 
for a parabolic problem is large, we can approach the time steps of an 
implicit scheme while retaining the advantages of an explicit scheme. 

Using an upwlnded scheme for a hyperbolic equation also results in 

eigenvalues that are not near the imaginary axis except for the longest of 

wave lengths. When one considers a pure one-sided scheme together with 

boundary conditions then the matrix is upper (or lower) diagonal and so all 
the eigenvalues lie along the negative real axis. Some experiments with this 
case are presented in [20]. 
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In our analysis we have assumed that one approaches a steady state fastest 
by using the largest possible time step. This is Intuitively obvious if one 
is using a time-like iteration procedure. However, when the iteration process 
is not time-accurate, the criterion for the optimal need not be connected 
with the time step. This is well known for A.D.I. type schemes. In section 6 
we shall show that when one uses residual smoothing, that the best strategy is 
not to choose the largest possible time step, even in one-space dimension. 
Recently, Jameson [9] has developed a multigrid code using the Runge-Kutta 
scheme as the smoothing operator. For this method one suspects that one would 
choose the so as to damp the high frequencies. Hence, one choice for the 
parameters , in one dimension, is to choose the a. so that 



for some 0 q, At and with the additional constraint that |G G| ^ 1, 

-IT ^0 ^ir. Since there is not a well-developed theory for the multigrid 
method for hyperbolic equations, it is not clear that (2.13) yields the 
optimal parameters. Computations indicate that one wishes some combination of 
large time steps together with good dampllng properties at high frequencies 
for the raultigrid procedure to be efficient. 
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3 . Time Step in Generalized Coordinates 

We consider the two-dimensional equation 

w + Aw + Bw = 0, (3*1) 

t X y 

where A and B are constant matrices which represent the gradient of the 
fluxes appearing in (2.1). We shall only consider Runge-Kutta methods in time 
and second-order central differences in space. The effect of the artificial 
viscosity on the time step is ignored (see previous section for a more 
complete discussion). From (2.5) we see that for any multistep Runge-Kutta 
method that the stability criterion is of the form 


At • d < K, 


(3.2) 


where d is the maximum (over 0 and spectral radius of 

D = A sin 0 + B sin (J). The constant K depends on the parameters of the 

Runge-Kutta scheme. 

The Euler equations in general coordinates (^>T)) are 

(JW)^ + = 0, (3.3) 

where 


p 


pq 


pr 

pu 

f = 

pqu + p 

8 = 

pru - P 

pv 


pqv - p 


prv + p 

E 


q(E+p) 


r(E+p) 


J = x_ y 
I Ti 


( 3 . 5 ) 
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Let 



-(Y-1)u 

pc 


-(y-1)v 

pc 


Y-1 

pc 



p 

0 


0 


1 

P 


0 


(3.12) 


0 




-(Y-1)u 


-(Y-1)v 


Y-1 


Then 

0 
0 
0 
w 

where a,b are given by (3.8) and w by (3.7). Hence, 

d = |w| + / a^ + b^ c. 


Dq = TDT 


.-1 _ 


w 

ac 

be 

0 


ac 

w 

0 

0 


be 

0 

w 

0 


(3.13) 


(3.14) 


Letting A? = Ati = 1, then the time restriction is of the form 


At < 


K*J ^ W 

+ /!?+? c |q| + |r| + / x^+y^4x^+y^ + 2 |x^x^+y^y^| c 


(3.15) 
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where r and q are given by (3.6). In the code FL052ST the stability 
criterion is implemently slightly differently as 


At 


K‘J 


|q| + |r| + (/ Xp + yp + / + y^)< 


(3.16) 


For most exterior problems a highly stretched mesh is used. In these 
cases the time step is governed by the area of the smallest cell which is much 
smaller than the area of the cells in the far field. To avoid this difficulty 
we use a different time step in each zone. This destroys the time accuracy of 
the solution but accelerates the convergence to the steady state. This local 
time step was first used in [13]. 

To see the effect of using a local time step we consider the radially 
symmetric wave equation 


u = — (ru ) . 
tt r r r 


(3.17) 


Discretizing the time variable we have 


lim 

At+0 


n+1 o n , n-1 

(At)^ 


7 

r r r 


(3.18) 


Let At be the local time step (^^)j rewrite (3.18) as 


lim 

At -*-0 


n+1 o n . n-1 « 

u. - 2u. + u. 2, V 

J 3 3 _ c (r) 




(ru ) , 
r r’ 


c(r) = 


(At), 




min 


(3.19) 


1 
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Therefore, using a local time step is equivalent to introducing an artificial 
wave speed that increases as one goes to the far field. Hence, the further a 
wave goes towards infinity, the faster it goes. As an example, we consider a 
mesh that increases exponentially with r. Then is proportional to r 
and so c(r) = r. Then (3.19) becomes 


“tt ■ 


(3.20) 


Let s = log(r), (3.20) becomes 


u. . = u . 
tt ss 


(3.21) 


If we only allow outgoing waves, then the solution to (3.21) is 


u = f(s-t) = f(log(r) - t). (3.22) 

Hence, if we begin with a wave of compact support at r = 1, t = 0, then at 
time tQ, the wave is centered at r = exp(tQ), i.e., the wave has moved 
exponentially fast. 


4. Artificial Viscosity 

The Runge-Kutta method with central differencing in space has two 
difficulties. The first difficulty is that the highest frequency is not 
damped, i.e., for a linear problem the neighboring points decouple. This odd- 
even decoupling prevents the possibility of driving the residual to machine 
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zero. With the nonlinear equations the variables are averaged at the cell 
faces before forming the fluxes. This nonlinearity couples all the neighbors 
together. Nevertheless this coupling is weak and convergence to a steady 
state can be slow. In order to accelerate the convergence a fourth-order 
linear viscosity is added to each equation. 

A second difficulty with central differences is that it does not enforce 
an entropy decrease across shocks. As such the central difference may 
converge to the wrong solution. We attempt to enforce the entropy condition 
by introducing an artificial viscosity. The fourth-order viscosity does not 
help near shocks. In fact there are theoretical indications that the fourth- 
order viscosity can be a destabilizing factor [14]. These observations are 
confirmed by numerical experiments. Hence the fourth-order viscosity is 
turned off in the neighborhood of shocks. In order to prevent oscillations at 
the shock an additional viscosity is added which behaves as a second 
derivative. This is in the spirit of the Navier-Stokes equations and seems to 
yield shocks without any overshoots. This viscosity is a nonlinear viscosity 
so that the formal order of accuracy is not affected by the addition of the 
artificial viscosity. 

We wish the artificial viscosity to accomplish two contradictory purposes. 
We want these terms to accelerate the convergence to a steady state. This 
implies that we should choose the viscosity as large as possible without 
violating any stability restrictions. On the other hand, we wish the 
viscosity to be as small as feasible so as not to affect the accuracy of the 
solution. Formally, the viscosity is of higher order than the truncation 
error. However, for a finite mesh, increasing the viscosity will decrease the 
accuracy of the steady state. 


T 
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In section 7 we will introduce an upwlnded scheme that does not require 
the use of an artificial viscosity. The use of an artificial viscosity has 
both advantages and disadvantages. The basic disadvantage is its 
artificiality. An artificial viscosity is not aesthetically appealing. 
Furthermore, there are invariably constants which must be adjusted for each 
case. This makes the code less robust. An upwind scheme has a built-in 
viscosity which should automatically adjust Itself to each situation. The 
advantage of an artificial viscosity is its flexibility. We can use the 
freedom of arbitrary constants or functions to tailor the code to an 
individual problem. This requires more work on the part of the user but it 
can be beneficial. The upwind schemes are more automatic. However, in cases 
where one does not want any viscosity, it is difficult to turn off since it is 
built into the algorithm. 

In order to introduce the artificial viscosity into the Runge-Kutta scheme 
we modify (2.2) and get 


where Lw represents the approximation to the Euler equations. is given 

by 


128 


) + D D_^ D w)l, 

^ +ri -n '• At +q -n ’ 


(4.2) 


and V 2 will be presented later. At is the local time step while J is the 
area of the zone. is an arbitrary constant (usually about 0.2) while 

and D_ represent forward and backward differences respectively. 5 and n 
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are the curvilinear coordinates and we have assumed A? = An = 1. As one 
approaches the boundary, can no longer be generated by (4.2). One 

alternative is to extrapolate the variables to an artificial position across 
the boundary and then to use (4.2). This is equivalent to using a one-sided 
approximation to (4.2). The viscosity added by is a completely 

dissipative mechanism. If we look at the differential level (4.2), (with 
V 2 = 0) leads to an approximation of 


w^ = Lw - (Kw ) , 

t XX XX 


K = 


C4 J 
Tt“ • 


(4.3) 


Multiplying (4.3) by w and integrating over all space we get 


Tit |w|^ dx = / (w, Lw)dx - / w(Kw^)^^ dx. (4.4) 


Integrating the last integral by parts twice and ignoring boundaries we get 


Y ^ / w^ dx = / (w, Lw)dx - / K(w^^)^ dx. 


(4.5) 


Hence, as long as K is positive the viscosity terms decrease the total 
energy. An alternative to (4.2) is to use 



+K 



w) + D_^ 

+n^At +n 



(4.6) 


This version contains both dissipative and dispersive characteristics. The 
same technique can be used on the finite difference level using summation-by- 
parts. One can then also Include boundary terms to see their effect. 
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Eriksson and Rizzi [4] choose boundary conditions so as to maximize the 
dissipation of the artificial dissipation. We shall see later that this may 
introduce errors into the steady state approximation. 

The form of V 2 is given by 

''a "]• 

0 is a switch which measures the gradients of the flow. In smooth regions 

2 

0 should be of the order of A while near shocks 0 should be of order 
unity. Hence, this viscosity is of fourth-order in smooth regions and does 
not affect the order of the scheme. 0^ is given by 


0 


X 

3+ V 2 


.k 




(4.8) 


Two possibilities for 9 are 



ic 


|Pj+l,k ~ Pj-l,k| 

|Pj+l,k ^Pj,k Pj-l,kl ’ 


K ~ 1.0, (4.9a) 


or 




IPj-H.k ~ ^Pj,k Pj-l.k' 


Ipj+iA ■ Pj.kl 




0.05. (4.9b) 


As with the fourth-order viscosity, (4.7) cannot be used directly at the 
points next to a boundary. Salas (see [1]) suggests extrapolating the 
variables to a virtual point outside the domain and then using (4.7). 
Eriksson and Rizzi [4] again Implement the boundary terms so as to Increase 
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the dissipative effect* However, computations indicate that the viscosity 
V 2 introduces errors near the boundaries* This error consists of two 
parts* Near the leading edge false entropy is generated which then forms an 
entropy layer along the body* At the same time a pressure jump occurs at the 
trailing edge which violates the Kutta condition* This problem is especially 
noticeable at high angles of attack* FL052ST was used to find the solution 
about a NACA 0012 at 10® angle of attack with a free stream Mach number of 
0*3* The solution is completely subsonic and so the Euler and potential 
solutions should agree* With the standard code a noticeable pressure jump is 
generated at the trailing edge* Since the flow is subsonic one can set V 2 = 
0* Without the second— order viscosity the pressure is continuous at the 
trailing edge* Furthermore, comparisons of the lift between the potential and 
Euler codes show that the Euler code underpredicts the lift by about 9% on a 
coarse 64x16 0 mesh On a finer mesh th^ lift is improved but is still much 
worse than the corresponding prediction of the potential code* Removing the 
V 2 component of the viscosity improves the lift prediction though it is still 
5% too low on the 64x16 mesh* Careful checks show that the pressure jump at 
the trailing edge is due to the tangential component of the V 2 viscosity. 
Hence, this difficulty is not caused by the difficulty in evaluating (4*7) 
near the boundary* The cause is that the switch (4*9) becomes large near the 
leading and trailing edge* Hence, these regions are treated as if there were 
a shock and a loss of entropy is created* This entropy layer extends over the 
length of the airfoil and is different on the upper and lower surfaces* 
Hence, the stagnation pressure on the upper and lower surfaces are different 
leading to a pressure jump at the trailing edge. If the tangential component 
of the V 2 viscosity is set equal to zero then the pressure jump disappears* 


T 
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However, there is no improvement in the lift prediction. Thus the lift is 
more sensitive to the normal component of the viscosity# 

In [3] it is suggested that the viscosity V 2 be multiplied by a linear 
factor which is zero near the body and one in the far field# This improves 
the accuracy in many cases, however, it is ad hoc# An alternative is to 
multiply the viscosity by M^# In stationary flow, shocks only occur in 
supersonic regions and so the viscosity is not changed near shocks# However, 
there is a stagnation point at the trailing edge and so the viscosity is 
turned off near the trailing edge# 


5# Enthalpy Damping 

With a second-order system in time one can add artificial terms that 
depend on the first-time derivative# Though such terms destroy the time 
accuracy they do not affect the steady state solution. These terms can be 
chosen so as to speed up the convergence to a steady state# Such applications 
have been used for SOR or the full potential equation# However, for a first- 
order system one cannot, in general, add on lower order terms since they are 
not zero in the steady state# For the Euler equations it is known that the 
total specific enthalpy is constant along each streamline in the steady state# 
If all the streamlines originate from a constant reservoir then the enthalpy 
will be constant in the entire region even in the presence of shocks# This 
constant enthalpy is known a priori from the inflow boundary condition# 
Therefore, one can add artificial terms to each equation that depend on the 
deviation of the local enthalpy from the steady state enthalpy# Such a 

forcing term is zero in the steady state and we will show that it can 
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accelerate the approach to the steady state. For simplicity of presentation 
we shall assume a one-dimensional Isentroplc fluid. We therefore consider the 
following modified equations. 


Pt + ^ 9 ^ + pu^ = <1) 

(5.1a) 

2 

u. + uu + p =0 

t X p ^x 

(5.1b) 

(|) = -ap(h - hp). 

(5.1c) 

2 2 

p Y - 1 ^ 2 » 

(5. Id) 


differentiating these equations, one obtains 


3^^ + 2up ^ + (u - c )p + u p + p^ u + 2uu p 

tt Xt XX t X t X X 


(5.2) 


- pu - p(— 1 P = <})^ + U<|) 
X '• p-'x X t 2 


also 


2 

^t + U(})^ = - a[(h - hQ + c )pj. + u(h - hQ>p^]. 


(5.3) 


We ignore all terms Involving the product of derivatives and then freeze the 
coefficients* We then have 


ptt 


2up 


xt 


- (c^ - 


u )p 


XX 


+ a(h - hQ + c )pj. - ocu(h - hQ)p^ 


0 . 


(5.4) 
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This is a convective wave equation in which a multiplies the first time 
derivative of p. This is similar to the acceleration procedure used for the 
full potential equation. To reduce this to a standard wave equation we 
introduce new Independent variables 

T) = X, T= + / 1 - t. M=-. ( 5 , 5 ) 

! 2 

c/ 1 - M 

where we assume that M £ 1. With this change, (5.4) becomes 

- (c^ - u^)p [d - 2M^(h - h ) + (1 - M^)c^]p^ 

/i-m2 (5.6) 

- au(h - h^)p^ = 0 
or 

Ptt - (5.7) 

Using a finite difference scheme we choose K proportional to d/Aq so as to 
reach a steady state rapidly. In terms of the original variables this implies 

aAt proportional to = ^ 5 =- . (5.8) 

(1 - 2M^)(h - hp) + c^(l - m"^) 

Note? Since o depends on 1/At the enthalpy damping is not a low order 

term and so it affects the stability limit. 

As noted above this procedure is valid only for subsonic flow. For the 
potential equation it is well known that adding a (|)^ term is not advisable 

in supersonic regions. Using enthalpy damping for the Euler equations it has 

been found experimentally that using a small amount (relatively) of enthalpy 
damping in the supersonic regions can accelerate convergence, especially for 
problems which are mainly supersonic. 
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2 

In the code FL052ST a = (1 - M ) + a 2 » a^, ot 2 constant. In the above 

derivation we used the primitive variables. For the conservative variables, 
(5.1) is replaced by 

= - ap(h - hp) 

(pu) j. + 1*2 ” “ 

(pv) j. + = - apv(h - hg), 


where represent the standard Euler terms. The finite difference 

approximation to (5.9) is 


- p“ = -AtL^ - aAtp“''’^(h" - h^) (5.10) 

or 


n+1 

P 


p” - AtLj^ 

' ■ • 
1 + aA(h" - hg) 


(5.11) 


We stress then in these formulas, a is always positive Independent of the 
sign of (h - hQ). Thus, the right-hand-side of (5.9) is not the equivalent 
of a forcing function that gives rise to an exponential decay in an ordinary 
differential equation. Replacing (h - hg) by |h - hQ| in the definition 
of <t> destroys the enthalpy damping and can lead to divergence. 

For the energy equation we append to (5.1) 

Sj. + = 0. (5.12) 


r~ ■ 
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Deriving the equation for E - p one gets a forcing term that depends on 
(h - hQ)^ which can lead to difficulties. In [8] this was fixed in a 

heuristic manner. An alternative is to replace (5.12) by 

h^ + = -8(h - hp) (5.13) 

A 

where 8 may depend on the Mach number. Let E = E - h^ p then combining 
(5.13) with (5.9) yields 

E^ = -(h - hQ)(aE + I p) = -[a(h - h^) + | ]e + | p. (5.14) 

We note that when 3 is large we force the enthalpy to be equal to its steady 
state value. Hence (5.14) is a generalization of the isoenergetic systems. 

Computations show that- the enthalpy damping is also useful in removing 
temporal oscillations. Using the enthalpy damping, with a Mach number 
dependence, yields a more monontone convergence to a steady state than in the 
absence of damping. This property is especially useful if one uses an 
acceleration procedure on top of the Runge-Kutta scheme. 


6. Residual Smoothing 

Residual smoothing was first Introduced by Lerat [12] for use with the 
Lax-Wendroff scheme. Jameson [9] later introduced a similar technique in 
conjunction with the Runge-Kutta schemes. We shall later compare the effects 
on both these methods. However, we first consider a two-step Runge-Kutta 
method. This scheme is given by 
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(1) n ^ ^ (n) 

u' = u + aAt Qu 


+ At Qu 


( 1 ) 


( 6 . 1 ) 


where Q denotes all the spatial derivatives. When an artificial viscosity 
is used it can be used either in both steps or else frozen at the same value 
for both stages. The residual smoothing is then given by 



where the product is over all space dimensions. One can also replace the 
operator on the left-hand-side of (6.2) with a full multidimensional elliptic 
operator. Since this operator only involves constant coefficients on a 
rectangular region (in computational space), it can be inverted by a fast 
solver. However, we shall see that one should not use large time-steps in 
conjunction with (6.2) even in one-space dimension. Since the difficulties 
are not concerned with splitting errors there is not much of an advantage to 
consider multidimensional operators that are not in split form. Numerical 
experiments indicate that a multidimensional Laplaclan is more effective in 
accelerating the flow to a steady state than (6.2) but not sufficiently fast 
to warrant the additional cost of inverting the multidimensional matrix even 
with a fast solver. 

We now consider the constant coefficient problem in one-space dimension. 
Let Q be the second-order central difference and Ignore the artificial 
viscosity. Taking the Fourier transform of (6.2) for u^ = u^^, we get 
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(l + 3 sin^ -l-jCG-l) = iX sin 6 - aX^ sin^ 0 
X = At/Ax 


(6.3) 


or 

_ _ , , iX sin 0 - aX^ sin^ 0 

t> - 1 + s— 5 . 

1 + 3 sin^ j 

The original two-step Runge-Kutta scheme (3=0) is stable when 


(6.4) 


X < 


/ 2a-l 
a 


(6.5) 


For general 3 , the scheme is stable if we choose 


a I /2 


3 > 2aX(X - . 


( 6 . 6 ) 


2 

Alternatively, if we replace 3 by aX , and a V 2 » o >. 2a, then the 
scheme is unconditionally stable. 

In three-space dimensions, (6.4) is replaced by 


(1 + aX^ sin^ |)(1 + oX^ sin^ |-)(l + aX^ sin^ f)(G-D = (6.7) 


with K = sin 0 + sin 0 + sin 5» For stability we require that G* G < 1. 
This occurs if and only if 
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-2a(l + a\^ sin^ sin^ |•)(l + aX^ sin^ -I) + X^ + 1 _< 0. 

( 6 . 8 ) 

A sufficient condition for stability is that 

-2a[ 1 + aX^(sin^ j + sin^ + sin^ !■)] + X^ + 1 < 0. (6.9) 

Hence, the three-dimensional scheme is unconditionally stable if 

- s a[sin 0 + sin ({> + sin 

o > max — — 5-5 ^ 5^* (6.10) 

p,<|),5 2(sin j + sin"^ ^ + sin^ j) 

Therefore, we have shown that a sufficient condition for the three-dimensional 
scheme (6.2) to be unconditionally stable is 

a V 2 > cr ^ 6a. (6.11) 

Hence, the three-dimensional version of (6.2) is unconditionally stable. 
Moreover, since the Runge-Kutta scheme gives a steady state which is 
Independent of At and (6.2) is in delta form we conclude that the steady 
state solution to (6.2) is also independent of At. We point out that the 
version of residual smoothing proposed by Lerat does not have a steady state 
independent of At. This is because the Thommen scheme used by Lerat has a 
At dependent steady state. Since the time steps used are not very large 
relative to the explicit time step this dependence on At may not be 
noticeable to graphical accuracyt; We also note that for the one-dimensional 
Lax-Wendroff two-step methods (Richtmyer, MacCormack, etc.) and for the two- 
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dlmenslonal Bursteln scheme that the solution after the Intermediate step Is 
independent of At. 

We wish to stress that even though (6.2) is unconditionally stable, 
choosing a large time step is not the best strategy. Even in one-space 
dimension choosing a very large time step severly retards the convergence to a 
steady state. This is not true, in one-space dimension, for the Lax-Wendroff 
scheme, if 6 depends on the matrix of the differential equation (see [12]). 
Hence, the use of residual smoothing with the Runge-Kutta scheme is 
fundamentally different than the backward Euler method. For the residual 
splitting method the inefficiency of large time steps has nothing to do with 
splitting errors. Let B = oX , (6.4) becomes 


G(6) = 1 + 


iX sin 0 - X^ sin^ 0 

1 ^ i ^ 

1 + oX sin Y 


a ^ V 2 
a > 2o 


( 6 . 12 ) 


For large X, (6.12) becomes 


G(0) 


2 2 
g sin 0 

® sin^ 0/2 


(6.13) 


We thus see that without artificial viscosity that the highest frequency, 

0 = IT is not damped. 

If we add an artificial fourth-order dissipation with coefficient v at 
each stage, then (6.12) is replaced by 


G(0) 


X(i sin 0 - V sin^ y) + X^(i sin 0 - v sin^ 

= 1 + ± i — 


1 + oX^ sln^ Y 


(6.14) 


As X gets larger, we get 
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2 

G(ir) ~ 1 + — > 1, (6.15) 

0 

and so the scheme is not stable for large X. We next consider the case that 
the same artificial viscosity is used at both stages. In that case (6.12) is 
replaced by 

iX sin 0[l + aX[i sin 0 - v sin^ —)] - Xv sin^ y 

G(0) = 1 + 5 ^ . (6.16) 

1 + aX^ sin^ Y 

As X approaches infinity, the coefficient of the artificial viscosity goes 
to zero relative to the denominator and so we recover (6.13) in the limit. 
Hence, when we freeze the artificial viscosity, the scheme remains 
unconditionally stable. However, when X is large, the highest frequency is 
damped only a small amount proportional to 1/X. Therefore, if we wish to 
minimize G, we do not wish to choose X large. An alternative is to xhoose 

a finite X so that / G (0)d0 is minimized. 

0 

We also note that if one adds the viscous terms from the Navier-Stokes 
equations, then a similar analysis holds. Hence, if one wants the scheme to 
be unconditionally stable then the Navier-Stokes terms should be frozen and 
one evaluated once per cycle. If the Navier-Stokes terms are reevaluated at 
each stage then the residual smoothing should be applied at every stage and 
not after the second stage. 

The above analysis was for a two-stage Runge-Kutta scheme. For a multi- 
stage scheme one should apply the residual smoothing after every even stage. 
If the total number of stages is odd then one should do an extra residual 
smoothing after the cycle is completed. One can also show that the residual 
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smoothing will not stabilize a one-stage Runge-Kutta with central differences 
for a hyperbolic problem but the method is unconditionally stable for a one- 
stage method for parabolic problems. In fact it gives the backward Euler 
methods. Hence, for the Navier-Stokes equations one should use residual 
smoothing after every stage. 


?• Highly Subsonic Flow 

It is well known that for very subsonic flow that the flow can be 
considered as incompressible. The use of the compressible fluid equations is 
considered as inefficient, since the fluid velocity is much smaller than the 
speed of sound. The use of an explicit scheme requires At to be bounded by 
1/c. However, the physical properties change over time scales of order 1/u 
which is much larger. Similar arguments hold for viscous flows with a high 
Reynolds number. Hence, it is usually agreed that explicit schemes are 
inefficient for highly subsonic flows. We shall now show that if one is only 
^®^®i^®sted in the steady state then a minor change to the code can greatly 
increase the efficiency of the explicit method. Even if one is using an 
implicit method the following changes should Increase the efficiency of the 
scheme since all waves have similar speeds. 

The Euler equations are expressed as 

»t ^x ®y “ 

where (x,y) represent general curvilinear coordinates (see (3.3)). Since we 
are only Interested in the steady state we replace (7.1) by the system 
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M^w^ + f +g =0. 
t X ®y 


(7.2) 


The requirements on M are that the matrix be nonsingular and that the 
original initial boundary value problem still be well-posed. It is 
straightforward to solve (7.2) with an explicit scheme. With an implicit 
method only the diagonal portion of the matrix to be inverted is changed. 
Though the code solves (7.2) we shall only analyze the constant coefficient 
problem 


M 


,-l 


w^ + Aw + Bw =0, 
t X y ’ 


(7.3) 


where the matrices M, A, B are constant and A = — , B = . 


,( 0 ) = 


„-l 


w'-' =Tw, A_=TAT, B„=TBT,M„ =TM T, where T 


,-l .-1 
0 


,.-l 


*0 - “ - » ''O 

(3.12). Then (7.3) can be converted to 


Let 

is given by 


w^®^ + A_ w^^^ + B- w^°^ = 0 
0 t Ox 0 y 


(7.4) 


with 


-• 



** 






q 

c 

0 

0 


r 

0 

c 

0 

c 

q 

0 

0 


0 

r 

0 

0 

0 

0 

q 

0 


c 

0 

r 

0 

0 

0 

0 

q 


0 

0 

0 

r 




.. 





*- 


(7.5) 


where q and r, defined in (3.6), are the contravariant velocity components. 


T 
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If Mq = 

I, then 

we have 

not changed the eigenvalues or the stability 

condition of (7.1). 

We 

o 

now consider the case that u 

+ v2 

' O 

' « c . We wish to 

-1 









choose Mq so that 

the eigenvalues of 

Mq Aq and 

Mq 

Bq are Independent 

of c. In addition, 

we wish 

Mq to be 

positive definite. This will imply 

that (7.3) is 

a symmetric hyperbolic system and so is 

well-posed. One 

choice 

is 



c^ 

0 

0 

0 















0 

1 

0 

0 

9 

(7.6) 



0 

0 

1 

0 






0 

0 

0 

1 












2 2 
where z = max(e , u 

+ 

v^). e 

is introduced so that 

the 

matrix 

is not 

singular at stagnation points 

• In particular, e 

= .Olc seems to give 

reasonable results • 

Transforming back 

to the curvilinear coordinates we 

define 

^ 9 











— u 

-V 

1 

<5^ 

u^ + v^ 



US^ 


-u2 

-uv 

u 

o 

2 


Q = 








(7.7) 


VS^ 


-uv 


V 

h 

c2 2 

= s2 



hs^ 


-uh 

"vh 

h 




Then 


m“^ = I + dQ, 


M = I + eQ, 


(7.8) 


where 


d = 1), 


^(4- O' 


(7.9) 
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We note that given the first row of Q, the following rows are derived by 
multiplying the first row by u,v,h, respectively. Hence, the product of Q 
times a vector requires only six multiplications. For use in a Runge-Kutta 
scheme we evaluate the fluxes in (7.2) as usual including the artificial 
viscosity. The vector of the changes in the variables Aw is then multiplied 
by the matrix M where the elements of Q are evaluated at the previous 

stage. The variables at the next stage can then be evaluated. 

9 2 
“X Z 

Let M - ~2 * largest eigenvalue of D = A sin 0 + B sin <j) is 

c 

given by 


X 


|w|(l + M ^) + J w^(l — ) + 4(a^ + b^)z^ 

2 


(7.10) 


where w, a, b are given in (3.8), (3.9). We see that at a stagnation point 
M = 0(e) and X = 0(e). For M = 1, X = |w| + / a^ + b^ c. Hence, at low 
Mach numbers the largest eigenvalues (and hence the time step) is independent 
of the sound speed c. At transonic sound speeds the largest eigenvalue is 
comparable to the case with M = I. Comparing with (3.16) we see that the 
preconditioned form (7.2) allows the use of a larger time step for all 
subsonic flows. Applications to the incompressible Navier-Stokes equations 
are presented in [19] together with extensions to supersonic flow. 


8. l^wlnd Schemes 

At each stage of the Runge-Kutta scheme we have replaced the flux 
derivatives by central differences. This necessitated the use of an 
artificial viscosity. This stabilized the scheme in both smooth regions and 


I 
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provlded an entropy condition in shocked regions. To avoid this artificial 
aspect one can replace the derivatives by appropriate upwind differences. 

Presently this is implemented using the flux splitting described by Van 
Leer [24] for the Euler equations. At each zone face 5,ri = constant, we 
calculate 

'^i+V2,j " '^i " (8.1a) 

w. (8.1b) 

6^ is constructed from the forward and backward differences using a switch. 
This switch prevents overshoots when the variables change dramatically over a 
zone width. The flux is then split into plus and minus contributions 
depending on the sign of the eigenvalues. These fluxes are calculated in a 
rotated system using the velocity components parallel and perpendicular to 
each coordinate direction. The result fluxes are then rotated back to yield 
the Cartesian fluxes. We then combine the plus and minus contributions from 
neighboring cells to obtain the flux at the face of the cell. This is done 
independently in each direction. Given the fluxes at all four faces of the 
cell we update the variables to the next stage. The coefficient of the Runge- 
Kutta scheme given by (2.7) were appropriate for central differences. Optimal 
values of the parameters for one-sided schemes are not known. Further details 
of the scheme are presented in [20]. 

Because of the extra logic involved in an upwind scheme it is more costly 
per stage than a central difference scheme. Though a central difference 
scheme requires an artificial viscosity it can be calculated once and then 
reused at each stage of the Runge-Kutta scheme. Since the viscosity of an 
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upwlnd scheme is built in, it is more difficult to perform some operations 
once and then reuse them during the k stage scheme. A four-stage scheme 
using upwind differences requires about three times as much computer time as 
the corresponding four-stage central difference scheme. In addition, the 
enthalpy damping technique of section 5 can not be used. This occurs since 
the enthalpy is not a constant in the steady state when using flux splitting. 
The stability limit for the one-sided scheme is only about two-thirds of that 
allowed for the central difference scheme. Hence, the present version of the 
upwind scheme is about five times slower than the central difference scheme to 
reach a steady state with a given tolerance. The chief advantage of the 
upwind scheme is its robustness. There is no need to choose constants for the 
artificial viscosity. Preliminary test indicate that the upwind scheme works 
for a larger range of inflow Mach numbers than the central difference scheme. 


9. Boui^dary Conditions 

In addition to advancing the scheme in the interior, it is necessary to 
supply boundary conditions. At the airfoil the normal component of the 
velocity is zero. Using the finite volume approach, it is only necessary to 
know the pressure on the airfoil. This pressure is found by the normal 
momentum equation. It is also necessary to give boundary conditions in the 
far field. When the flow is subsonic at infinity one should specify three 
conditions at Inflow and one condition at outflow. 

In one-space dimension, in order to decrease the energy as fast as 
possible and reach a steady state rapidly, one should specify characteristic 
conditions. Diagonalizing the wave equation one gets 


T 
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u 


t 


Vj. = V^, 0 < X < 1 

(9.1) 


u(l,t) = av(l,t), v(0,t) = Pu(0,t) 


u,v are the characteristic variables. Specifying characteristic boundary 
conditions is equivalent to a = P = 0. From (9.1) it follows that 

j- j (u^ + v^)dx = (a^ - l)v^(l,t) + - l)u^(0,t). (9.2) 

Hence, choosing a = p = 0 minimizes the right-hand-side. 

Therefore at inflow, one should specify three characteristic variables. 
However, we require that the enthalpy be constant over the entire field” in the 
steady state. To achieve this it is necessry to specify enthalpy at Inflow. 
This condition replaces one of the characteristic boundary conditions. For 
nonlinear problems it is more appropriate to specify the Riemann variables 
rather than characteristic variables. Hence, at Inflow we specify 

p/p^ 

V 

h=^. 

p 

u is the component of velocity perpendicular to the boundary while v is the 
component parallel to the boundary. The fourth boundary condition is found by 
extrapolation from the interior. For stability it is preferable to 
extrapolate characteristic variables [5]. Hence, at the inflow we extrapolate 


(9.3a) 

(9.3b) 

(9.3c) 



At the outflow boundary we reverse the procedure and specify 


and extrapolate 


2c 

(9.4a) 

p/P^ 

(9.4b) 

V 

(9.4c) 

h 

(9.4d) 


An alternative to (9.4) is to use nonreflecting boundary conditions [2], 

Numerical experiments indicate that the far field boundary exerts a much 
greater influence on the drag and lift coefficients than does the boundary 
condition at the airfoil. 


10. Conclusions 

We have discussed many of the components of the code FL052ST. This code 
gives a rapid solution to the Euler equations in both two- and three-space 
dimensions for a variety of geometries and a range of Mach numbers. 

With central differences it is easy to increase the accuracy of the 
differences to fourth-order or spectral accuracy. Using an 0 mesh all 
variables are periodic around the airfoil. Therefore, high-order differences 
do not encounter any boundary difficulties. Using spectral methods a Fourier 
scheme would be appropriate. With a C mesh all boundaries are in the far 
field. Hence, one can simply reduce the order of accuracy near the outer 
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boundarles where the flow is smooth enough for second-order accuracy to be 
adequate. The approximation normal to the boundary requires more care with 
regard to boundary conditions at the airfoil surface. Hence, it is reasonable 
to consider second-order differences normal to the airfoil while using higher- 
order methods parallel to the body. In this case one must be careful in 
approximating the metric derivatives so that the constant flow in the far 
field remains a solution to the finite difference equations. Work is also in 
progress on extending the code to the Navier-Stokes equations. 
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